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ABSTRACT 


The objective of this program was to develop a computerized mathematical model 
of the combustion response function of composite solid propellants, with particular 
attention to the contributions of the solid phase heterogeneity. Such a modGl would 
be useful to design propellant formulations that would minimize the tendency to drive 
combustion instability in given applications. Although it is established that 
anmonium perchlorate particle size distribution has a significant effect upon this 
driving, previous models have treated all or portions of the combustion zone as 
homogeneous . 

The one-dimensional model treats the solid phase as alternating layers of AP 
and binder, with an exothermic melt layer at the surface. Solution of the Fourier 
heat equation in the solid provides temperature and heat flux distributions with 
space and time. The problem is solved by conserving the heat flux at the surface 
from that produced by a suitable model of the gas phase. An approximation of the 
BDP flame model is utilized to represent the gas phase. By the use of several 
reasonable assumptions, it is found that a significant portion of the problem can 
be solved in closed form. A method is presented by which the model can be applied 
to tetramodal particle size distributions. 

A computerized steady-state version of the model was completed, which served 
to validate the various approximations and lay a foundation for the combustion 
response modeling. The combustion response modeling was completed in a form which 
does not require an iterative solution, and some preliminary results were acquired. 

If is concluded that the solid phase heterogeneity does per se influence the 
time lag and phase shift mechanisms responsible for combustion driving, and thereby 
the response function. Although the model was not fully evaluated by comparisons 
with experimental response function data, this conclusion is supported by the nature 
of the preliminary results. However, these results also indicate that the role of 
AP cannot be attributed to the solid phase alone — at least insofar as the solid 
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phase Is represented by this model. Potential deficiencies are Identified, and 
areas of future work are recommended. 
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SECTION 1 


INTRODUCTION 


Experimental data have established that ammonium perchlorate (AP) particle 
size has a significant effect upon the pressure-coupled response function of 
composite solid propellants (1-14). Moreover, the effect cannot be attributed 
simply to changes in burning rate or formulation; the effect appears to involve 
the composite propellant heterogeneity as well (15). Classical theories of combustion 
driving (16) have assumed a homogeneous propellant, and are therefore inadequate to 
explain the combustion instability characteristics of composite propellants. The 
cwimunity has come to rely upon experimental measurement of the combustion response 
in T-burners, and work in recent years has been devoted largely to improving the 
method (17-19). Although experimental measurement serves several important purposes, 
interest in the theoretical work has revived because the acquisition and interpre- 
tation of full complements of data continue to be expensive and do not furnish a 
phenomenological mechanism for the guidance of propellant R&D. 

Viewing the ccwnbustion zone as the region between the thermal wave penetration 
in the depth of the solid and the location of the flames in ttie gas, there are 
several ways in which the composite propellant heterogeneity can manifest itself. 

Two schools of thought have arisen: one which emphasizes the solid phase, treating 

the gas as a homogeneous source of propellant heating; and one which emphasizes the 
gas, continuing to treat the solid as a homogeneous medium. 

The solid phase proponents may be represented by Lengelle & Williams (20), 

Kumar (-21) and Cohen (13). Lengelle and Williams performed a one-dimensional analysis 
of a solid having sinusoidal thermal properties. Although the model was too idealized, 
it made the essential point that the heterogeneity augments the combustion response 
depending upon the periodicity of the thermal properties (and, therefore, particle 
size and spacing). Kumar introduced an exothermic surface melt layer, purpor«.c:!iy 
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representative of the AP surface, in an otherwise homogeneous solid. The most 
significant result of this model was a mechanism by which zero-exponent propel- 
lants could exhibit a positive combustion response and pressure effects.^ Since 

Kumar did not treat in-depth heterogeneity, particle size effects appeared in 
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the solid only through the effect of burn rate on melt layer thickness. Cohen 
postulated two characteristic parameters for the solid phase, each dependent 
upon particle size. One was a measure of frequency response, the other of 
thermal response. He assumed that all response function curves could be 
determined by the pressure exponent plus these two parameters, and used experi- 
mental data for a standard propellant as calibration. The essential result is 
that decreasing particle size increases the peak response and shifts the peak 
to higher frequency. He further assumed that multifnodel propellants could be 
treated by linear superposition of the constituent results, and predicted multi - 
peaked response functions. Although the method is not founded upon a formal 
analysis, it is being used to guide propellant tailoring {15, 22). Extensive 
applications have revealed some qualitative and quantitative deficiencies, for 
example a tendency to over-emphasize the effect of fine sizes (22). 

The gas phase proponents may be represented by Hamann (23), Beckstead (24), 
and Click and Condon (25). All utilize some form of the "BDP" model of steady- 
state combustion (26) to represent the gas phase details, and none consider the 
solid to be heterogeneous. Hamann performed a perturbation analysis upon the 
entire BDP model, but did not report any results. Beckstead used the BDP model 
to calculate values for the parameters v/hich are called for by the homogeneous 
theory of Denison & Baum (27). This approach of combining unrelated models is 

1 The classical the^ories produce a response function proportional to pressure 
exponent, yet it is well known that plateau (zero exponent) and mesa (negative 
exponent) propellants have exhibited combustion instability. The AP melt may 
be analogous to the foam zone of such propellants. 

2 An effect of burn rate on melt layer thickness would also appear in double- 
base propellants, so cannot be a sole basis for the role of AP particle size 
in composite propellants. Furthermore, AP size effects persist at constant 
burning rate when catalysts are used or solids loading or distribution are 
varied. 
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open to question (15), and a review of the technique reveals several deficiencies.^ 
Click & Condon employed a similar approach, but used a modified BDP model (28) 
for polydisperse particle size distributions and the method of Zeldovich & 
Novozhilov (29) as an alternative to the theory of Denison & Baum. Comparison 
of results with experimental data was disappointing. Considerable improvement 
was noted, however, when the Cohen postulates were incorporated into the method 
to position the peak response and peak resoonse frequency (25). 

The following consensus appears to emerge fr«n this background. First, there 
is a need to account for the melt layer and the in-depth solid phase heterogeneity 
of composite propellants. Second, there is a need to provide an analytical basis 
to test, confirm or modify the Cohen postulates. Third, the representation of the 
gas phase also should address the heterogeneity of composite propellants by em- 
bodying some form of the BDP model rather than the homogeneity of the classical 
theories. Accordingly, it was the objective of this program to develop an 
analytical model of the combustion response of composite solid propellants with 
particular attention to these contributions of the propellant heterogeneity. 


3 For example, calculation of the "A" and "B" parameters of the homogeneous 
theory by this method reveals much too small a pressure effect to account 
for measured changes in the response function curve with pressure. As 
with Kumar, particle size effects appear only through changes in steady- 
state burning rate properties. 
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SECTION 2 


MODEL PREMISES AND ASSUMPTIONS 

2.1 REPRESENTATION OF THE SOLID PHASE 

The solid phase is represented by a "sideways sandwich", following the 
concept of Lengelle & Williams, as shown in Figure 1. This picture, really, 
does nothing more than state that the analysis is one-dimensional, so its dis- 
similarity to real propellants is of no greater concern than is the use of a 
one-dimensional treatment. Such a treatment assumes that the lateral processes 
are negligible in comparison to the normal processes. The solid is considered 
to be semi -infinite, having alternating layers of AP and binder. The thickness 
of the AP layers is nominally equal to the particle size, with exceptions to 
be noted later. The thickness of the binder layers is equal to the interstitial 
spacing as determined from the statistical georetry (26, 30). The surface AP 
layer contains a thin melt layer, which is justified experimentally (31-33), 
and follows the model of Kumar as a region of exothermic reactions in accordance 
with an Arrhenius law. "^he melt layer is "thin" in that the melting point of 
AP approximates the surface temperature during deflagration (32, 33). The 
analysis is linearized for small harmonic pressure perturbations, and is con- 
cerned viith pressure-coupling only. 

Although the model contains the convective heating term to represent the 
regression of the material at the mean rate f, the geometry of the layers is 
taken to be fixed with the AP layer always at the surface. This assumption is 
similar to the statistically fixed geometry used in the BDP model. One conse- 
quence of this assumption is that the "pulsation" mechanism associated with moving 
layers is excluded. The pulsation mechanism is open to question in view of the 
aggregate macroscopic properties of the propellant, such that there is no 
coherence over the propellant surface for resonance, although Lengelle & Williams 
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Figure 1. Physical Model Representation 






offered an argument in its favor. In any event, the present physical model will 
not be conducive to such resonance because the thermal properties are not 
sinusoidal in space. Thus, the present model examines only the heterogeneity 
in relation to the thermal wave, which is a uniform property of the propellant 
and therefore considered more realistic at this point than pulsation. On the 
other hand, another consequence of this assumption is that in-d'*pth heterogeneity 
will not be important in those cases where the particle size is larger than the 
thermal wave (generally, coarse particles and high burn rate). Whether or not 
this consequence is unduly restrictive remains to be seen, but the fact that 
catalyzed coarse propellants are more stable than fine propellants would seem to 
permit it. The statistically fixed geometry requires the surface AP layer 
(including the melt) to have a thickness particle size ( 26 ). 
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2.2 REPRESENTATION OF THE GAS PHASE 

The function of the gas phase in the analytical scheme is to transform the 
oscillating pressure into an oscillating heat flux boundary condition at the 
surface of the AP melt. An approximate font of the BDP model nas been selected 
to represent the gas phase, and the gas phase is assumed to be quasi-steady 
(viz, the gas phase responds instantly to changes in pressure). The approxima- 
tion assumes a single flame above the propellant surface, where all gas phase 
reactions occur, as illustrated by Fig. 1. 

Presuming that the condensed phase heterogeneity has the dominair effect 
upon the response function of composite propellants, it would not seem to matter 
what particular model v/ere chosen to represent the gas so long as it provides a 
reasonable boundary condition. Thus, the Denison & Baum model could have been 
chosen to preserve some systematic order to the analytical development. This 
was not done for two reasons: First, it would have presumptuously ignored 

those developments addressing the gas phase heterogeneity in a BOP model frame- 
work. Second, the fact that the Denison & Baum model is heavily dependent upon 
fluctuations in flame temperature is coming to be viewed as a serious deficiency 
of that model. Variations in flame temperature are fourth order with respect 
to variations in pressure. On the other hand, variations in flame standoff , 
which are not addressed by Denison & Baum but which are a significant aspect of 
the BOP mods! (and a key to the p?,i*ticle size effects in the gas phase), are 
first order with respect to pressure variations. 

The approach, then, is a perturbation of an approximate form of the BDP 
model with respect to flame standoff as well as flame temperature. Perturbation 
of the model itself is considered to be proper, whereas use of the model to 
calculate parameters for substitution into a different model is open to question. 
With the approximate model, the task is much simpler thar that undertaken by 
Hamann . 
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2.3 REPRESENTATION OF MULTIMODAL PROPELLANTS 

Multimodal propellents are represented by adjacent columns of layers, 
as shown by Figure 2. Each particle i,»7e tops a column, corrected by the factor 

. The remaining AP layers consist of the finest size in the distribution. 
This assumption is justified by the fact that the finest size will fill the 
smallest interstices between coarser particles, and will be valid so long as 
the thermal wave does not penetrate to a subsequent coarse sublayer. The 
assumption oF limited thermal wave pe’-'^tration will be valid for practical 
propellants because of the influence of the fine size in raising the burning 
rate. The binder layers are all of the same thickness, as computed by the method 
of Ref. (30). All of the columns are constrained to the same mean burning rate, 
which implies one effective flame height or characteristic size for the dif- 
fusion flame. However, the columns are allowed to have different response 
functions and the aggregate value for the propellant is a weighted average from 
the constituent columns. 

2.4 APPROACH 

The analysis is performed in two parts. First, the model is derived in its 
steady-state version. The steady-state version serves to calculate mean values 
required for the time-dependent model, verify that the boundary conditions will 
be reasonably accurate as measured by the ability to reproduce experimental (or 
formal BDP) results, estaLnjh the zero-frequency (no oscillation) limit, verify 
som of the model assumptions, and provide an initial check on the method of 
solving for the termperature profile in the solid. The steady-state model is 
one of the subroutines of the computerized response function model. Secondly, 
the time-dependent model for the combustion response function is derived. In 
calculating the response function, experimental values of mean burning rate are 
used as input in order to minimize the effects of uncertainties in the st-^ady- 
state modeling. 
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Figure 2. Layered Model of the Solid Phase for a Trimodal Propellant: 

Sizes Dj, D^. D^. Shaded Layers are Binder. The Top of the 

Surface AP Layer is a Thin Melt 


SECTION 3 


THE STEADY-STATE MODEL 


3.1 SOLID PHASE EQUATIONS 

The heat conduction equation in the melt layer is written as: 

•'a f = «a‘>a‘*a*®>''> 

where k = thermal conductivity 
p = density 
c - heat capacity 
W = weight fraction 
Q = heat of decomposition 
subscript a = denotes AP 

A = kinetics prefactor 
E = activation energy 
R = gas constant 
r = burning rate 
T = temperature 
X = distance into the solid 
The boundary conditions are: 

* -''a S = ^sS'-'Tw - ^ 

* ' ’'m- -"a ^ ' ^o> * «a»s^n, 

where subscript w = denotes wall or surface 
subscript s = denotes mean propellant 
subscript m = denotes melt 
subscript o = denotes the deep solid 


(1 


(2 

(2 
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Equation (1) may be written in dimensionless form as follows: 


d^T . dx _ Q 

^* 37 - B exp 


-E/Rt 




where t = (T-To)/T^-Tj,) 
y =Ix 

K = thermal diffusivity of the medium, k, for Eq, (3) 

Q 


B = A ^ ri, ^ 
r 


w2 "a ca 


H = 


WQ 


In the thin melt layer, x*^ 1, so Eq. (3) may be approximated as: 


+ ^ = C exp [-D(l-x)] 


where C 
D 


Bexp (-E/RT^) 


RT. 


2 <Tw'o> 


w 


The boundary conditions become: 


surface: y = 0 , ^ = -Z, ( 1 +Hj + H^) 


melt/crystal: y - y^^^, dx _ 


-Z, {r + 
a m m 


where = 


p$Cs 


^ ^a"a 


Making the following transformation: 


r - 4- 


n = D(l-x) 
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Eq. (4) becomes: 


(c-l-n/D) ^ exp (-n) 


( 8 ) 


r 

Recognizing r:/D «) where and letting K = sr, Eq. (s) becomes: 




(9) 


Eqs. (5a & b) become: 


n=0. C = =1-ZJ1+H^ + HJ 


’w a 


(10a) 


n=n„» ^ = T,„ “ Z, (t +H„) 

m m m a m m 


(lOb) 


Eq. (9) may be integrated in closed form. Applying Eqs. (10) and some algebra 
yields an expression for burning rate in the following form: 


K = - 


(V-0- 


1 - exp 


( 11 ) 


Making all substitutions, burning rate is expressed as a function of wall 
temperature and other constants as follows: 

2 


2Wa 


RT 


w 


-2, C^rTM-To)E(Tw-To) 

> y I I «. 


exp(-^ ) l-exp( 


LU_ 


-rT* 

X‘w Iw 


^ - j(^)(l-Z,) - Z, ^ 


( 12 ) 


- 1 


This relation is to be matched with a relation between burn rate and wall 
temperature from the gas phase model. 

Analysis of the steady-state problem and, as will appear later, soli’ .ion of 
the time-dependent problem, also require a description of the steady-state 
temperature profile in the solid. The temperature profile in the melt layer is 
determined from Eq. (9). Integration and application of Eq.(ioa) yields: 


C = 1 ± + 2 [Kexp(-n)“Cf^] 


(13) 
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where C„ = 

The physics of the temperature decay requires the negative root. Substituting 
the definitions for c and t yields: 

^ = -n + dVi+ 2 [Kexp(.ii)-C„] (14) 

An order of magnitude analysis shows that the square root term will always 
be much larger than n in the melt layer. Thus, Eq. ( 14 ) reduces to the form: 


dn 


“^a + bexp(-n) 


= Ddy 


(15) 


where a = l-2(^ 


b = 2K 

This equation can be integrated in closed form for y(n); it is transcendental 
as n(y). 


y 


1 

D^a 


In 


Hz -Vz^ 2z 
"-Jd2+2d 


(16) 


where d = 2a/b 
z = dexp(n) 

The melt layer thickness may be calculated by evaluating Eq. (16) at n=n^. 


For all layers beneath the melt layer, the right-hand-side of Eq. (3) vanishes. 
The gradient condition at the top of each layer is similar to Eq. (5b), without the 
heat of fusion. 


y 


V — = 
^Top* dy 



( 17 ) 


\»here Z * Z_ when entering an AP layer 

Q 

Z » Zj^ when entering a binder layer 
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Thus the temperatures in each layer beneath the melt layer follow the recurring 
form: 

’"'Top ■ ^ ’Top (18) 

where Ay = thickness of the particular layer, and uses 

the appropriate k. 

For a homogeneous propellant, Z=l, so Eq. (18) properly reduces to the resuH for 
a homogeneous propellant. Note also that the temperature and the gradient will 
properly tend to zero together, as y approaches infinity. 


3.2 GAS PHASE EQUATIONS 

If it is assumed that all reactions occur at the flaire height, the heat 
conduction equation has the form of Eq. (1) with a zero right-hand-side. Taking 
the flame to be at x=0 and the wall at the flame height (x=x*), the temperature 
distribution is: 


T - Tf-(T^-T„) 


l-exp(-c) 

l-exp(-c*) 


where c = ux/kq 

u = gas velocity normal to the surface 
subscript g = denotes gas 
subscript f = denotes flame 


(19) 


For convenience, c is set equal to y by employing the continuity relation, 
PgU=p$r, and assuming that the ratio of heat capacity to thermal conductivity 
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IS equal for gas and solid . Then the gradient at the wall may be written, from 
differentiation of Eq. (19), as: 


y=y' 


dy 


t w' I -exp 


where y 


>*- 


rx* 


( 20 ) 


4. This assumption is comparable to the general assumption that these properties 
are temperature-insensitive. The conductivity of propellant gases may be 
calculated, but is not well-known. 
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At the same time, the gradient at the wall must satisfy the energy requirements 
of the condensed phase: 


where subscript b = denotes binder 

Eq. (21), when compared to Eq. (2a), states that the binder heat of decom- 
position is being positioned on the gas side of the wall. The binder does not 
appear at the surface in the model of Fig, 1, but does in reality exist over 
portions of the surface. The model surface is the AP melt. Therefore, in the 
framework of this model, the decomposing binder is external to the melt layer 
so may be represented by a heat absorption on the gas side. Eq. (21) may be 
re-written as: 

y = y*. =(T„-v + F (22) 

where F = 


Equating Eq. (20) and Eq, (22) yields: 

Tw = (Tf + F - Tq) exp (-y*) - (F-T^) 

which is the required matching relation in burn rate and wall temperature. 

The remaining unknown, x*, is determined by an approximate fit of the 
effective flame height from the BDP model: 



where D-j = particle size 
p = pressure 

C- = 24.6 for r in cm/sec, p in atmospheres, D,j 
in microns and x* in microns. 


(23) 


(24) 
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SECTION 4 


THE TIME-DEPENDENT MODEL: CALCULATION OF THE RESPONSE FUNCTION 
4.1 SOLID PHASE EQUATIONS 

Using dimensionless quantities as previously defined, the time-dependent 
heat conduction equation in the thin melt layer is written as: 


Eq. (25) is the time-dependent form of Eq. (4). Denoting mean values as 
barred and perturbed values as primed, and employing Eq. (4) to describe the mean 
portion, Eq. (25) may be written as: 

“ Ifet) ^ (-n^[exp(DT')-l] (26) 

3y 

where tu = frequency of oscillations 
n - 10 

For harmonic perturbations, t' = exp{iwt), the time-dependent term of 
Eq. (26) may be re-written to provide an ordinary differential equation: 


2 

^ - iiJx' = Cexp(-n) [exp(Dx')-l] (27) 

dy"^ 


If x' is of second order and p sufficiently small, Eq. (27) may be approximated 


as: 



CDt' exp(-S) 


(28) 


18 



The problem is linearized to small perturbations, and is therefore rt?stricted 
to linear, instability. Variations in burning rate are of the order of variations 
In pressure, but variations in surface temperature are second or third order with 
respect to variations in burning rate. Thus, for second order pressure pertur- 
bations, t' is at most of third order. The exponential in n is approximately 

2 

unity in the melt layer, and the product CD is of the order 10 . This would 
permit frequencies as high as lOKHz to satisfy the approximation for cases of 
interest. Since the quasi-steaoy assumption for the gas phase model will also 
restrict the frequencies to values less than lOKHz, use of Eq. (28) is satisfied 
by that assumption. 

Applying the boundary conditions at the surface must recognize that the 
surface is fluctuating relative to the mean surface (y=0) position. Since the 
fluctuating burning rate will also be of the form exp(ia>t) in the linearized 
problem, it follows that the surface position is given by: 


>S = 


ini 


(29) 


The boundary conditions are: 


'S’ 


t' = T. 


W 


W ‘ "b> 

r 


(30a) 

(30b) 


where g' = the gradient in t’ at the surface, on 

W 

the negative (gas) side of the boundary. 

Eq.. (30b) is derived by substituting Eq. (2a) into Eq. (21), and perturbing 
In the non-dimensional form. For small perturbations, the quantities at the 
actual surface may be related to the quantities at the mean surface as: 
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y = yg. 


(31a) 


t' » t' 


y=0 


1 r* - 




^•v=0 


dT’ . n.l 

W ' ^ 'y=o 



(31b) 


where g = the gradient in t 


Eq. (28) is integrated numerically. The functions exp(-n) and g are available 

in terms of y from Eq. (16). The integration provides the mean surface values, 

which are then converted to the actual surface values by means of Eqs. (31). The 

actual surface values must satisfy Eqs. (30); g' is determined from the gas phase 

w 

model, and t' is related to r' by perturbing Eq. (12). The perturbation of Eq. (l2) 

W 

yields: 

- Vjt; (32) 


where Vg = (1+ |-) x f 
0 = E/(RT^) 

X 


l-2x(l-T m) 

exp[D(l-T^)]-l 


1 


Although the solution appears to require iteration, it will be shown sub- 
sequently that it does not (cf. Subsection 4.3). 

For all layers beneath the melt layer, the right-hand-side of Eq. (27) 
vanishes so the time-derivative term cannot be neglected. Thus: 




iftx' = 0 


(33) 
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The boundary conditions at each interface are generally; 


y=y 


rop, 



(34a) 


dx* „ ^a.b , 

dy Stop 


(34b) 


where the first Z subscript is used when entering an AP layer and the second 
is used when entering a binder layer. Assuming that the product pc for the 
melt is equal to that for the solid AP, the ratio of Z does not appear for 
that first layer of solid AP which joins the melt. The g'j^p refers to that 
value of the perturbed gradient on the negative (upper) side of the boundary, 
which drives the behavior below. The however, is preserved on the 

positive side. Eq. (34b) is similar in form to Eq. (30b), but there is no phase 
change heat beneath the melt layer and the sublayers do not oscillate relative 
to the mean surface. Eq. (33) may be integrated to produce the following 
recurring formulas: 

X x' - a' 

Top Z. «Top 

t'st'jqp exp(A^Ay) * [ exp(X 2 Ay)-exD(x^Ay)] (3j) 

Xg-Xi 


dx' _ 
dy 




^1^‘top 


Ilzk 

f 

~b,a 


g' 


Top 




[X2exp(X2Ay)-\^ exp(x^Ay)] (3 


where x^ * * f ‘ ^ V^+4in 


X 2 “ “ ^ ^ ^l+4in 

Each layer uses the appropraite ratio of Z, and the appropriate k for ay, x^ 
and X2. 
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Eqs. (35) and (36) properly reduce to the results for a homogeneous solid 
(Z=l , and the difference between g' and vanishes throughout leaving only 
the firrt term on the right-hand-side). However, these equations retain the 
growing exponential term (in Xg, of positive real part) for the layered solid 
aecause the formulation of the problem does not Impose explicitly the boundary 
condition at Infinit/ upon each layer. This condition is the requirement 
that T* and di'/dy vanish together. For the homogeneous solid, this condition 
would inmediately set the exponential in equal to zero. For the layered 
solid, it can only be said that a layer will eventually be reached where this 
condition may be approximated for all practical purposes. Since the effect of 
the heterogeneity on the perturbations also disappears when the perturbations 
disappear, this approximate condition may be expressed in the form of the 
homogeneous solid; 

t' -► 0, X^t' (37) 

In other words, there is a depth below which the propellant can be treated as 
homogeneous for purposes of the perturbation problem. Consider that this occurs 

XL. 

below the N— AP layer. Then Eq. (37) may be expressed as: 


y=y» 


T “€ 


(38a) 


dT* , 
dy *^ls ^ 


(38b) 


where c = an extremely small number 

x^j » value of x^ using for the homogeneous region. 
Eq. (38b) represents the perturbed gradient at the top of the homogeneous 
region,, on the positive (homogeneous) side of the boundary with the N~ AP layer. 
This may be converted to the perturbed gradient on the negative (AP) side of the 
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boundary through use of Eq. (34b), recognizing that Z=1 In the honiDgeneous 
region. Thus, Eqs. (38/ become: 


y 



(39aj 


^TcpN “ h ^Is*^ 

Where subscript TcpN ' 

layer, which would be the "top" 
condition for the succeeding homogeneous 
region were the calculation to continue. 

Thus, the problem of the solid phase Is now properly closed. Eqs. (3Q) 
define the conditions at the "top", and Eqs. (39) define the conditions at 
the "bottom". The conditions at the bottom are now known, but the conditions 
at the top depend upon the gas phase. Additional discussion with respect to the 
solution for the solid phase is deferred to Subsection 4.3. 


4.2 GAS PHASE EQUATIONS 

The perturbed form of Eq. (24) is simply: 



Since the analysis will take into account flame temperature perturbations, 
a relation between the flame temperature perturbations and the wall temperature 
perturbations Is required. T'.is is obtained by perturbation of Eq. (23). 

Usli.^ dimensionless quantities as previously defined, Eq. (40) for the flame 
standoff perturbations, and linearizing the exponential In perturbed quantities, 
there results: 

exp(y*) - (1+Tp)y*exp(r'')[?- - 2 f- ] (41) 

P r 
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where xp 


F 




0 


The perturbed gradient on the gas side of the wall Is derived from an 
energy balance. The energy being transmitted Into the solid is the difference 
between the heat release In the gas and the energy required to raise the gas 
temperature to the flame temperature: 

->^9 ^ - S'VV 

where Qp = heat release in the gas {negative for an exotherm). 

Using the same relation between gas and solid thermal properties as led 
to Eq. (20), Eq. (42) may be written in dimensionless form as follows; 


wh6*"6 g = gradient in x at the wall on the negative (gas) side 

W 




Of 

(Tw-^o) 


The perturbed dimensionless form of Eq. (42) Is: 



— [H, + 1)] + (t.' 

r ’ 



) 


Hp may be eliminated by combining the relat 


vrts In Eqs. (43), (5a) and (21)^. 


5 For this purpose, Eq. (21) may be written as: dx/dy = Z(g + H. ). This 
Is the steady-state analog of Eq. (30b) , and comes from substituting 
Eq. (2a) into Eq. (21) In dimensionless form. 


Substitution into Eq. (44) yields: 


9„' = - ^ (Hip) 

Substitution of Eq. (41) into Eq. (46) yields, after combining terms: 

g ' = — (1+Tp)(2y*exp(y*)-1) + r ‘ (exp(y*)-l ) > ■Ml+xr)y*exp(y*) 
w - r w r 


Eq. (47) is the necessary matching relation for Eq. (30b). Thus, the 
formulation of the time-dependent problem is complete. 

Eq. (45) also may be substituted into Eq. (43) to yield: 

% " -n+Tp) 

Eq. (48) is the dimensionless form of Eq. (22), so consistency is verified. 

4.3 SOLUTION FOR THE RESPONSE FUNCTION 

By combination of Eqs. (30b), (31a), (31b), (32) and (47), it is possible 
to derive an explicit t ression for the response function in terms of solid 
phase constants, gas phase constants =»nd one key parameter characteristic of 
the solution for the solid phase. This key parameter is the ratio of the 
perturbed graident to the pe*turbed temperature at the mean surface: 


The components of this ratio appear individually in Eqs. (31), and the 
ratio is defined here as K^. In the classical homogeneous theory^ this ratio 
is always . In this work, the ratio will depend upon the modeled heterogeneity 
and A 2 as well as on A-|. 

The expression for the response function is: 



Where R = response function 

V 3 = (see Eq. [5a]) 

Vg = (see Eq. [32]) 

C = (see Eq. [4]) 

VgA = (exp(y*)-l) (see Eq. [47]), coefficient of t^') 

VgB = (1+Tp)y*exp(y*) (see Eq. [47] , in coefficients of 
P 7 p and r'/r). 

If K 2 depends upon g^', the problem will require an iterative solution. 
However, it turns out that K 2 is an intrinsic property of the solid phase, 
independent of the surface boundary condition, just as is X.j for the homogeneous 
solid. The reason is that the solid phase is described by linear homogeneous 
differential equations. All solutions of such equations of second order are 
of the form: 

f = C^f,(X) + C 2 f 2 (X) 
where f = denotes functions 

C-| = constant associated with surface boundary condition 
C 2 = constant associated with deep solid condition 
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Thi;- form does not take into consideration the relative importance of the com- 
ponents of the BDP multiple flame structure and their differing deoendencies 
upon the parameters of Eq. (24). However, it is a reasonable representation 
selected for mathematical convenience and with the perturbation analysis in mind. 

After substitution of Eq. (24) into Eq. (23), Eqs. (12) and (23) are two 
equations for the unknowns r and T . These equations are solved by iteration. 

W 

Once r and are known, it is possible to calculate the melt layer thickness 
and the thermal profile in the solid. 

The steady-state model is not intended to be used to calculate burning rates 
for multimodal propellants. That would require some definition of an effective 
particle size or flame height for use of Eq. (24). Rather, experimental values 
of burning rate are used to determine the wall temperature from the solid phase 
model, whence the effective flame height and particle size may be determined from 
the gas phase model. This is the general method for determining the mean values 
for use in the response function model. An option is provided to predict these 
values for unimrdal porpellants, and is used in this work only to validate as- 
pects of the modeling. Those results are discussed in Section 5. 

Note that the effects of particle size appear only as gas phase effects in 
the steady-state model. This conforms with a generally-accepted view of steady- 
state combustion (34). Although Eq. (12) (from the solid phase) may influence 
the magnitude of the particle size effect, particle size does not appear in that 
equation. The thermal profile in the solid does not enter into the calculation 
for mean burning rate. However, this does not preclude the importance of solid 
phase effects in determining the role of particle size in the time-dependent 
model for the response function. 
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The gradient is: 

g = C^f^'(x) + C 2 f 2 '(x) 

It is required that g/f go to in the deep solid. This provides a 
relation for C 2 : 

f^'(x)-x^f^(x) 

Substituting this expression for Cg into the ratio g/f, it is verified 
that C-| cancels out. Of course, g/f will vary from case to case and will 
vary with x for a given case; the important point is that it does not depend 
upon C-j . This property was observed numerically in the course of computer 
program development, which was based initially upon an iterative scheme. 
Accordingly, it suffices to go through Eqs. (35) and (36) in the layers, 
beginning with Eqs. (39)®, and then solve numerically through the melt layer, 
just one time for . Knowing Eq. (50) for R is solved by complex 

o 

arithmetic . 

It is of interest to examine the form of Eq. (50) in comparison to the 
form obtained from the homgeneous theory (Ref. 16). The latter can be 
written as: 


_ nnAB 

fi[AB-(l+A)] +iA-X^[fi-iA] 

where n = pressure exponent 

A = a solid phase parameter, not to be confused with 
A as defined in Eq. (1 ) 

B = a gas phase parameter, net to be confused with 
B as defined in Eq. (3) 

X-| = as defined in Eqs. (35, 36); it is of opposite 
sign in Ref. (16). 


6. Subroutine LAYRSP of the computer program 

7. Subroutine MLTLRP of the computer program. 

8. Subroutine GLIFP of the computer program. 
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It Is noted that Eqs. (50) and (51) are identical in form and, to some 
extent, they are similar in substance. The parameter Vgg combines condensed 
phase and gas phase terms, as does (pAB) of Eq. (51). The parameters C and 
Vg, "1 though different, are condensed phase terms as is A of Eq. (51). The 
parameter Vg^ is a gas phase term, as is B of Eq. (51). The parameter Vg is 
related to C, and therefore, a part of the analogy to A; but important 
differences from A derive from the finite melt layer. The parameter K 2 is 
the heterogeneous analog of x^. The parameter is purely a consequence of 
the heterogeneity, so there would be no analog for it in Eq. (51). The 
identity in form, and the similarity in substance, suggests that the hetero- 
genity as described by this model will not produce radical changes in the 
qualitative behavior of the response function. 

The ze»*o-frequency limit of Eq. (51) is the pressure-exponent, n. 

Pressure exponent does not appear explicitly as such in the steady-state 
model described herein. Nevertheless, it is of interest to examine the zero- 
frequency limit of Eq. (50). It is readily apparent that a non-zero response 
function at zero frequency requires that the following relationship be satisfied: 


Lim Ko = “at 


The satisfaction of Eq. (52) is verified by combining Eqs. (31), (32), (5a), 
the derivative of Eq. (15) applied at y=0, and a perturbation of Eq. (2a). In 
general, must be solved numerically, but at zero frequency it ■$ possible 
to derive an expression which reduces to -C/V^. Since Eq. (52) is satisfied, 
the indeterminate form of Eq. (50) that results may be evaluated to yield the 
non-zero response function: 


Lim R 

Q 


'6B 


2''6B 


(53) 
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Eq. (53) may be thought of as an "effective" pressure exponent, extracting an 
implicit property of the model. Note that it depends upon condensed phase 
terms as well as gas phase terms. Numerically, it is found to be consistent 
with pressure-dependence as calculated from results of the steady-state model. 
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SECTION 5 


MODEL RESULTS 


5.1 THE STEADY-STATE MODEL 

Results oV the steady-state model have been obtained in order to evaluate 
some of the important model premises. The essential results are tabt'^ated in 
Table I. 

Table I presents various results, compared with data and with BDP model 
results, for A-13 propellant used as a standard case. The first set of results 
compares burning rate as a function of pressure. The model results compare 
very well with the data. It should be emphasized, however, that these results 
should not be construed to imply that this model is “better than" the BDP 
model. The second set of results compares surface temperature. Experimental 
values of surface temperature are reportedly in the neighborhood of 850°K (26). 
It is observed that this model produces higher surface temperatures than the 
BDP model, and a somewhat greater sensitivity to pressure. However, the results 
are reasonable. The results also confirm the assumption that variations in 
surface temperature are second order (or smaller) with respect to variations 
in burning rate. The third set of results compare flame standoff distance. 

This model uses one flame. The BDP results are for the primary flame (sum 
of diffusion and reaction heights), and the values to the right of the slash 
are for the AP monopropellant flame when that flame moves closer to the surface 
than the primary flame. When that happens, the BDP model employs an energy 
partitioning which may be thought of as some single flame having an effective 
height between the two shown. On that basis, the flame heights from this model 
are roughly a factor of 3 greater than from the BDP model but the qualitative 
behavior with pressure is the same. The value of Cp in Eq. (24) was adjusted 
to achieve good agreement with the burning rate data; values of other constants 
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TABLE 1: COMPARISCWS OF BURNING RATE AND OTHER 

QUANTITIES FOR A-13 PROPELLANT 
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500 19.9/18.3 35.7 

700 19.5/5.7 30.2 

900 19. A/3.1 25.8 

1100 19.3/3.1 23.6 


are the same as used in the BDP model. It is concluded that, for purposes of 
the time-dependent analysis, the model conforms reasonably well with steady- 
state reality. 

The ability to reproduce measured effects of AP particle size on burninq 
rate was tested with a series of propellants analogous to A-13. These propel- 
lants were the subject of a low pressure L* instability study oerformed by 
Ramohalli (35). The comparison of burning rates at 100 psia is as follows: 


PARTICLE SIZE(u) DATA (cm/sec) MODEL (cm/sec) 


40 

0.41 

0.38 

90 

0.27 

0.28 

200 

0.23 

0.19 

360 

0.19 

0.13 


Again, the agreement is reasonable. 

There are two other aspects of the steady-state model results which merit 
discussion: the melt layer and the heterogeneity in relation to the thermal 
wave. 

The melt layer thickness is computed to be of the order of microns or 
less, which is consistent with experimental observation and the thin melt layer 
assumption. Its dependence upon heating rate involves a tradeoff between sur- 
face temperature and the steepness of the thermr.l gradient. Theoretically, it 
will disappear at such low burning rate that the surface temperature does not 
reach the AP melting point, and also will approach zero at very high burning 
rate where the gradient is very steep. Althcugh the layer is thin, it was 
considered improper to neglect it for mathematical convenience because the 
characteristic time of its dimension corresponds to high frequencies of 
interest. 

For particle sizes in excess of 40y, the thermal wave will not ot-netrate 
the first AP layer under conditions of interest. The implication is that, ex- 
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cept for the melt layer, the solid can be considered homogeneous in determining 
its role. However, this is not true for the fine sizes which are generally 
utilized in practical propellants. An estimate for a 2\i AP propellant reveals 
that, at 1000 psi, the temperature doesnotfall to within 10% of the bulk 
temperature until about 5 pairs of AP -binder layers are traversed. Further, 
if the 2\i AP is a component of a multimodal propellant, the burn rate will be 
lower such that the thermal wave will penetrate more layers of the column con- 
sisting of the 2m AP. As a result, it appears that the role of solid phase 
heterogeneity will be limited to melt layer heterogeneity in the intermediate- 
coarse size regime, but that in-depth heterogeneity can be important in the 
fine size regime. This distinction is one consequence of the present fixed- 
geometry model; were the layers permitted to move to evoke the pulsation 
mechanism, then the in-depth heterogeneity would always be important. The 
distinction was considered signifcant in view of the experimental importance 
of fine AP (15). 
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5.2 THE TIME-DEPENDENT MODEL 

The effects of the solid phase heterogeneities are most likely to appear at 

combinations of fine AP and low burning rate. Therefore, a test case consisting of 

a Z\i AP propellant at a burning rate of 0.47 cm/sec was selected for evaluation. 

Except for the particle size, this last case would correspond to A-13 propellant 

at 300 psi. Results are shown in Fig. 3. The solid line is for the heterogeneous 

propellant. The long-dash line is for AP and binder thermal properties equal to 

mean propellant thermal properties; therefore, it is for melt layer heterogeneity 

g 

only, the propellant below the melt layer being homogeneous . The short-dasn 
line is for a completely homogeneous solid; i.e., the distributed heat release in 
the finite melt lyacr is now concentrated at the surface only and the entire region 
beneath the surface is homogeneous^^. It is observed that the effects of the 
heterogeneities are small. 

Figure 4 compares results for a 90y AP propellant, which is A-13 propellant, 

with the Figure 3 results. Thus, the effect of particle size at a constant 

burning rate is shown. In the framework of this model, a constant burning rate 

implies a constant wall temperature and constant dimensionless flame properties; 

thus, any difference is due to solid phase heterogeneities. An effect of the 

heterogeneity does appear, but again, it is small. It is noted that the results 

for A-13 are virtually identical to the homogeneous solution displayed in Fig. 3. 

In the case of A-13, the thermal wave does not penetrate the surface AP layer 

and the melt layer thickness is about 1% of the particle size; thus, the solid is 

homogeneous for all practical purposes. In the case of the 2y propellant, the 

thermal wave penetrates 15 AP layers and the melt layer thickness is about 1/3 

of the particle size; thus, the solid is heterogeneous, but the effect of the 
9- Z-l. 

10- Vj = C = Z = 1 ; Kg = ; Vg = e-x- This would correspond to the Denison and 

Baum model except for differences in the modeling of the gas phase, and 
differences in the values of combustion constants due to the use of this model 
(including the finite melt layer) to reproduce steady-state burning rates. 
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Figure 3. Response Function Calculations for a 2 fiAP Analog of A- 13 Propellant 





heterogeneity appears to be small. Significantly, the effect is small with respect 
to peak response frequency as well as magnitude. Thus, it can be concluded that 
the expected effects of the heterogeneity are not being represented by this model. 

The theoretical results for the A-13 propellant, moreover, do not agree with 
experimental data. The A-13 is a JANNAF standard propellant, and has been well- 
character izad in T-ourner testing {Ref. 36). The experimental peak response is 
about 4 and occurs at rough\/ 300 Hz. Thereafter, the response declines to a 
value of approximately 1.5 at 1000 Hz. These theoretical “esults show a peak 
response of 0.58 at 125 Hz, and a value of 0.51 at lOOOHz. Therefore, the theory 
shows a relatively slight peak and at too low a frequency. Seen theoretical re- 
sults are a consequence of the combustion parameters, analogous to the "A" and 
"B" (or ''a'') parameters of the Denison and Baum model. Presumably, the agre^ement 
with data could be improved by selection cf a different set of parameters. However, 
a ground rule of this study was that the parameters would be predetermind b. 
considerations of a credible -teady-state model. Given that mcdel, it is in- 
appropriate to change the ccnstoPts arbitrarily. According to Ref. (27), a 
relatively large value of the ‘'\i" 'or “a") paramete' ‘'ill constrain the time-dependent 
model to the prediction of small peaks, and a relatively small value of the "A" 
parameter will produce low peak response frequencies. Sy analogy, that is the 
situation here. Since it would be inappropriate to juggle parameters, it must be 
concluded that there is a mechanistic deficiency in the time-deoendent model. 

Figure 5 presents theoretical results showing the effect of burning rate 
for a constant particle size. The 2y AP propellant was selected for this illust'^a- 
tion. The higher burning rate is representative of this propellant; the lower 
burning rate may be thought of as a suppression for purposes of Figure.s 3 and 4. 

The effect on peak response frepi ;-ncy demonstrates further that the modeled 
heterogeneities are of little consequence. It is observed that the peak response 
frequency varies nearly with the square of the burning rate, which is the result 
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1.23 Io./mc 



Figure 5. Effect of Burring Rate at Constant AP Particle Size for the Zft Analog of 
A- 13 Propellant 



for homogeneous propellants. On the other hand, the Fig. 3 and 4 results show that 
the effect of the heterogeneity on peak response frequency is about 10%. According 
to the Cohen postulates for the effect of heterogeneity, the peak response frequency 
in Fig. 5 should vary with the first power of burning rate, and in Fig. 4 it should 
vary inversely with particle size. These effects are not being produced by this 
model. 

Figure 6 compares theoretical results with experimental data for a bimodal 

propellant. Although the shape of the theoretical curve is reasonable, the peak 

response magnitude and freauency are again underpredicted. Also, the zero frequency 

% 

limit is overpredicted, reflecting a deviation from the measured steady-state 
pressure exponent. It appears that further work is necessary in order to implement 
a propel mechanism for the combustion response. 
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Figure 6. Comparison of Theory and Experiment for SP-540 Propellant at 500 psi 





SECTION f 

CONCLUSIONS AND RECOMMENDATIONS 


An analytical model has been developed which incorporates mechanisms of 
solid phase and gas phas^ heterogeneities into the calculation of steady-state 
and linear time-dependent combustion properties of composite solid propellants. 
Although the model satisfactorily describes the steady-state combustion properties, 
it is deficient in describing the time-dependent combustion response characteristics 
in several respects. Use of a consistent set of combustion constants produces 
peak response magnitu(?^s and frequencies which are too low in comparison to 
experimental data, auJ which are not significantly affected by the AP particle 
size per se. Altnough an effect of the solid phase heterogeneities is predicted 
by this model, the evfect is so small quantitatively as can be neglected in future 
work. Therefc^LS the role of AP cannot be attributed to the solid phase alone 
unless some c mechanism is incorporated into the theory. It is recommended 
that the concept of moving layers be re-examined, including justification for the 
coherence of such a mechanism. It is further reconmended that the perturbed BDP 
model be examined to represent the heterogeneity of the gas phase. It is desired 
not only to achieve the effects of the heterogeneity, but also to justify a set 
of values of the combustion constants that will properly position the response 
function curve. It appears necessary to modify both the solid phase and gas phase 
models in order to achieve those purposes in a consistent manner. 
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